Community Sampling

In [87]:

    

from collections import namedtuple
from itertools import cycle

from os.path import join, exists
from os import makedirs, chdir,getcwd,fchmod
from glob import glob
import shutil

import time

import numpy as np

import pandas as pd

import re

import scipy as sc
from scipy.stats import norm
from scipy.stats import expon
import scipy.interpolate as interpolate
import scipy.integrate
from scipy.fftpack import fft, ifft

from lmfit.models import Model
from lmfit import conf_interval

%matplotlib inline
import matplotlib
import matplotlib.pyplot as plt

from Bio import SeqIO
from Bio.SeqRecord import SeqRecord

from Bio.Seq import MutableSeq
from Bio.Alphabet import IUPAC

import docker

import configparser

import seaborn as sns
sns.set_context("paper")


    

In [44]:

    

def create_mount_run(image,mount_dir,cmd,envs):   
    container = cli.create_container(
        image=image, command=cmd, volumes=['/mnt/vol'],
        host_config=cli.create_host_config(binds={
            mount_dir: {
                'bind': '/mnt/vol',
                'mode': 'rw',
            }
        }),
        environment=envs
    )
    ctr=container.get('Id')
    cli.start(ctr)
    cli.wait(ctr,60*60*24*10)
    return cli.logs(ctr)


    

In [45]:

    

def Gekv(B,C):
    return 2**(1 - B)/(C*np.log(2))*(1 - 2**(-C))


    

In [46]:

    

def Gekv(B,C):
    return 2**(1 - B)/(C*np.log(2))*(1 - 2**(-C))


    

In [47]:

    

class Community:
    "A community defined by genome references (Biopython SeqRecords) and corresponding growth parameters."
    def __init__(self,name,acc,growth_param,td,mapper,env):
        self.name = name
        self.conf = local_conf(join(td,name),mapper)
        self.d_conf = local_conf(join('/mnt/vol',name),mapper)
        self.args0=['/home/Programs/PTR-Pipeline/ptr_pipeline.py','-c',join(self.d_conf['node_path'],'project.conf'),'-e','hedani@chalmers.se']
        self.env=env
        
        self.create_dirs(True)
        save_config(self.conf,self.d_conf)
        self.fetch_ref(acc)
        
        self.pop = self.init_population(acc,growth_param)
        self.distribution = self.community_distribution()
        self.samples = []
        
    def init_population(self,acc,growth_param):
        records=open_records(glob(join(self.conf['ref_path'],'Fasta','*.fasta')))
        population = namedtuple("population", "B C D l seq cells")
        
        #refId=[x.id for x in records]
        #acc=keys(comm)
        #refInd=[acc.index(x) for x in refId]
        
        pop = {}
        
        #for i,rec in enumerate(records):
        # add [:-2] to a if . removed 
        for i,a in enumerate(acc): 
            pop[a]=population(B = growth_param[i][0],C = growth_param[i][1],
                             D = growth_param[i][2], l = len(records[a]),
                             seq = records[a], cells=growth_param[i][3])
        return pop
    
    def ptr(self):
        for i,a in enumerate(self.pop.keys()):
            print a+": "+str(2**growth_param[i][1])
    
    def community_distribution(self):
        d=np.array([Gekv(p.B,p.C)*p.l*p.cells for p in self.pop.values()])
        return d/d.sum()

    #def ab_(self):   
    #    return 
    
    def sample(self,nr_samples):
        nr_samples=np.array(nr_samples*self.distribution)
        nr_samples=nr_samples.astype(np.int)
        
        samp_run=[]
        for i,p in enumerate(self.pop.keys()):
            samp_run.append(inverse_transform_sampling(self.pop[p].C,nr_samples[i],self.pop[p].l))
        
        self.samples.append(samp_run)
    
    def write_reads(self):
        for i,samp in enumerate(self.samples):
            write_reads(samp,self,self.conf['data_path'],self.name+str(i))
    
    def compare_fit(self):
        if not self.samples:
            print "The community is not sampled, please run community.sample(nr_samples)"
            return;
        
        err_hfit=[]
        err_pfit=[]
        res_fit=[]
        
        for i,samp in enumerate(self.samples):
            for acc in self.pop.keys():
                try:
                    depth_file=join(self.conf['output_path'],self.name+str(i),'npy',acc+'.depth.npy')
                    best_file=join(self.conf['output_path'],self.name+str(i),'npy',acc+'.depth.best.npy')

                    signal=2**(np.load(depth_file))
                    signal=signal/(signal.sum()/len(signal))#*self.pop[acc].l)
                    from_ptr=np.load(best_file)

                    res=fit_signal(signal,self.pop[acc].l)
                    res_fit.append(res)
                    
                    err_hfit.append((self.pop[acc].C-res.best_values['C'])/self.pop[acc].C)

                    err_pfit.append((self.pop[acc].C-from_ptr[2]+from_ptr[3])/self.pop[acc].C)
                    print "Simulated value: "+str(self.pop[acc].C)
                    print "Error from this fit: "+str(err_hfit[-1])+ ' value: ' + str(res.best_values['C'])
                    print "Error from initial PTR fit "+str(err_pfit[-1])+' value: ' + str(from_ptr[2]-from_ptr[3])
                except Exception as Ex:
                    print "Ex"
                    pass
        return [res_fit,err_hfit,err_pfit]
    
    
    def fetch_ref(self,acc=''):
        acc_path = join(self.conf['node_path'],"acc")
        f = open(acc_path, "w")
        
        if not acc:
            f.write("\n".join(self.pop.keys()))
        else:
            f.write("\n".join(acc))
        f.close()
        
        create_mount_run('centos/hedani',td,self.args0+['fetch-references','-s',join(self.d_conf['node_path'],'acc')],self.env)

    def build_index(self):
        create_mount_run('centos/hedani',td,self.args0+['build-index'],self.env)

    def run_pipeline(self):
        create_mount_run('centos/hedani',td," ".join(['/bin/bash -c "']+self.args0+['make']+[';']+self.args0+['run"']),self.env)
        
    def collect(self):
        create_mount_run('centos/hedani',td,self.args0+['collect'],self.env)
        
    def create_dirs(self,from_init=False):
        if exists(self.conf['node_path']) and from_init:
            shutil.rmtree(self.conf['node_path'])
        
        for d in [self.conf['node_path'],self.conf['data_path'],self.conf['output_path'],self.conf['ref_path']]:
            if not exists(d):
                makedirs(d)

    #def complexity(self):
        #len(self.pop)
        #return complexity


    

In [48]:

    

import numpy as np

def cartesian(arrays, out=None):
    """
    Generate a cartesian product of input arrays.

    Parameters
    ----------
    arrays : list of array-like
        1-D arrays to form the cartesian product of.
    out : ndarray
        Array to place the cartesian product in.

    Returns
    -------
    out : ndarray
        2-D array of shape (M, len(arrays)) containing cartesian products
        formed of input arrays.

    Examples
    --------
    >>> cartesian(([1, 2, 3], [4, 5], [6, 7]))
    array([[1, 4, 6],
           [1, 4, 7],
           [1, 5, 6],
           [1, 5, 7],
           [2, 4, 6],
           [2, 4, 7],
           [2, 5, 6],
           [2, 5, 7],
           [3, 4, 6],
           [3, 4, 7],
           [3, 5, 6],
           [3, 5, 7]])

    """

    arrays = [np.asarray(x) for x in arrays]
    dtype = arrays[0].dtype

    n = np.prod([x.size for x in arrays])
    if out is None:
        out = np.zeros([n, len(arrays)], dtype=dtype)

    m = n / arrays[0].size
    out[:,0] = np.repeat(arrays[0], m)
    if arrays[1:]:
        cartesian(arrays[1:], out=out[0:m,1:])
        for j in xrange(1, arrays[0].size):
            out[j*m:(j+1)*m,1:] = out[0:m,1:]
    return out


    

In [49]:

    

def pxn(x, C, l):
    return (2**(1 + C - (2*C*x)/float(l))*C*np.log(2))/(float(l)*(-1 + 2**C));


    

In [50]:

    

def pxn2(x, C):
    return (2**(1 + C - (2*C*x))*C*np.log(2))/((-1 + 2**C));


    

In [51]:

    

def piecewise_prob(x,C,l):
    conds=[(0 < x) & (x <= float(l)/2), x > float(l)/2]
    funcs=[lambda x: pxn(x,C,l)/float(2),lambda x: pxn((l-x),C,l)/float(2)]
    return np.piecewise(x,conds,funcs)


    

In [52]:

    

def piecewise_prob2(x,C):
    conds=[(0 < x) & (x <= float(1)/2), x > float(1)/2]
    funcs=[lambda x: pxn(x,C,1)/float(2),lambda x: pxn2((1-x),C,1)/float(2)]
    return np.piecewise(x,conds,funcs)


    

In [53]:

    

def inverse_transform_sampling(C,n_samples,l=1):
    x = np.linspace(0, float(1)/1.7, 100)
    y = np.linspace(pxn(.5,C,1),pxn(0,C,1),100)

    pdf=interpolate.interp1d(x,pxn(x,C,1))
    cdf = [scipy.integrate.quad(lambda x: pdf(x),0,i)[0] for i in np.linspace(0,.5,100)]
    inv_cdf = interpolate.interp1d(cdf, np.linspace(0,.5,100))
    
    r = np.random.rand(n_samples)
    v=l*np.round(inv_cdf(r[:int(n_samples/2)]),4)
    v2=l*(1-np.round(inv_cdf(r[int(n_samples/2):]),4))
    v=np.concatenate([v,v2])
    
    return v.astype(np.int)


    

In [54]:

    

def read_pair(header,ind,seq,pos,l):
    def circular_yield(x,pos,sub,l):
        if (pos>l):
            return x[pos-l:pos-l+sub]
        elif (pos+sub>l):
            r=pos+sub-l
            x2=x[pos:pos+sub-r]
            return x2+x[0:r]
        else:
            return x[pos:pos+sub]

    #base_error_rate = .02
    #mutation_rate = .001
    #nuc=["A","C","G","T"]
    
    # generate a normally distributued insert size of mean 300 and sd 40 
    #insert_size = int(norm.rvs(size=1,loc=300,scale=30)[0])
    insert_size = int(500)
    
    r1 = circular_yield(seq,pos,100,l)
    r2 = circular_yield(seq,pos+insert_size,100,l)
    
    # flip reads according to seqs error and mutation rate
    #ra=[]
    #for r in [r1,r2]:
    #    rr=np.random.random(100)
    #    m_ind=[i for i,j in enumerate(rr) if j < base_error_rate]
    #    char=np.random.choice(nuc,len(m_ind))
    #    a=list(r)
    #    for i,ind in enumerate(m_ind):
    #        a[ind]=char[i] 
    #    ra.append("".join(a))
        #r_tmp=r_tmp[:m_ind[0]]
        #for i,index in enumerate(m_ind[1:]):
        #    r_tmp = r_tmp[:index] + char[i] + r[index + 1:]
    #[r1,r2]=ra

    #nrs=[str(np.random.randint(low=1000,high=2000)),str(np.random.randint(low=10**4,high=2*10**4)),str(np.random.randint(low=10**5,high=2*10**5))]
    rec1=SeqRecord(r1,id=header+"_"+str(pos)+"."+str(ind),description="RANDXXLAB"+str(pos)+"_"+str(pos+insert_size+100)+":0:0"+"/1")
    rec2=SeqRecord(r2,id=header+"_"+str(pos)+"."+str(ind),description="RANDXXLAB"+str(pos)+"_"+str(pos+insert_size+100)+":0:0"+"/2")
    #rec1=SeqRecord(r1,id=header+"."+str(ind),description="FCB00YLABXX:6:"+nrs[0]+":"+nrs[1]+":"+nrs[2]+"/1")
    #rec2=SeqRecord(r2,id=header+"."+str(ind),description="FCB00YLABXX:6:"+nrs[0]+":"+nrs[1]+":"+nrs[2]+"/2")
    
    #rec1=SeqRecord(r1,id=header+"_"+str(pos)+"_"+str(pos+insert_size)+"/1",description="")
    #rec2=SeqRecord(r2,id=header+"_"+str(pos)+"_"+str(pos+insert_size)+"/2",description="")
    rec2=rec2.reverse_complement(id=header+"_"+str(pos)+"."+str(ind),description="RANDXXLAB"+str(pos)+"_"+str(pos+insert_size+100)+":0:0"+"/2")
    #rec2=SeqRecord(r2,id=header+"."+str(ind)+"/2",description="")
    rec1.letter_annotations['phred_quality']=[17]*100
    rec2.letter_annotations['phred_quality']=[17]*100
    
    #rec1.description=
    
    return(rec1,rec2)


    

In [55]:

    

def open_records(fasta):
    records={};
    for f in fasta:
        handle = open(f, "rU")
        tmp=list(SeqIO.parse(handle, "fasta"))
        m=re.match('.*(N[CTZ]_([A-Z]{2})*[0-9]{6}).*',tmp[0].id)
        records[m.group(0)]=tmp[0].seq
        handle.close()
    
    return records


    

In [56]:

    

def write_reads(samples,comm,directory,name):
    f1 = open(join(directory,name+"_1.fastq"), "w")
    f2 = open(join(directory,name+"_2.fastq"), "w")
    
    for i,p in enumerate(comm.pop.keys()):
        for j,pos in enumerate( samples[i].tolist()):
            r1,r2 = read_pair(name+"_"+p,str(j+1),comm.pop[p].seq,pos,comm.pop[p].l)

            SeqIO.write(r1, f1, "fastq")
            SeqIO.write(r2, f2, "fastq")
    
    f1.close()
    f2.close()


    

In [57]:

    

def plot_samples(samples):
    y=[np.bincount(s) for s in samples]
    for p in y:
        plt.plot(p/float(p.sum()),marker='.',linestyle='None')


    

In [58]:

    

def local_conf(td,mapper):   
    return {
            'project': 'ptr_simulation',
            'cluster': '',
            'job_name': 'ptr_simulation',
            'job_nodes': '1',
            'cpu_cores': '4',
            'estimated_time': '',

            'node_path': td,
            'ref_path': join(td,'References'),
            'data_path': join(td,'Data'),
            'output_path': join(td,'Out'),
            'doric_path': join(td,'DoriC'),

            'mapper': mapper,
            'ref_name': 'sim',
            'nr_samples': '1',
            'samples_per_node': '1',
            'email': 'hedani@chalmers.se',
            'data_prefix': '',
            'start_ind': '1',
            'job_range': '1',

            'ftp_url': '',
            'data_url': ''
        }


    

In [59]:

    

def run_community_analysis(exec_string,td,mapper,acc,growth_param,community,nr_samples,topname):
    chdir(td)
    comm=[]
    docker_td='/mnt/vol'
    
    for i,c in enumerate(community):
        accs=[acc[x] for x in c]
        grs=[growth_param[x] for x in c]

        wd='comm_'+str(topname)+str(i)
        docker_fd=join(docker_td,wd)
        fd=join(td,wd)
        
        # create dirs
        if 'd' in exec_string:
            if exists(fd):
                shutil.rmtree(fd)
            for d in [wd,join(wd,'Data'),join(wd,'Out')]:
                if not exists(d):
                    makedirs(d)
        
        chdir(wd)
        lconf=local_conf(fd,mapper)
        conf=local_conf(docker_fd,mapper)
        save_config(lconf,conf)
        
        args0=['/home/Programs/PTR-Pipeline/ptr_pipeline.py','-c',join(conf['node_path'],'project.conf'),'-e','hedani@chalmers.se']
        
        if 'f' in exec_string:
            f = open("acc", "w")
            f.write("\n".join(accs))
            f.close()

            create_mount_run('centos/hedani',td,args0+['fetch-references','-s',join(conf['node_path'],'acc')])

        if 's' in exec_string:
            # create community and sample
            # name,acc,growth_param,td,mapper
            comm.append(Community(wd,accs,grs,td,'bowtie2'))
            comm[i].sample(nr_samples)
            comm[i].write_reads()

        if 'b' in exec_string: create_mount_run('centos/hedani',conf['node_path'],args0+['build-index'])

        if 'm' in exec_string: create_mount_run('centos/hedani',conf['node_path'],args0+['make'])

        if 'r' in exec_string: create_mount_run('centos/hedani',conf['node_path'],args0+['run'])

        #chdir('..')
        
    return comm


    

In [60]:

    

def fit_signal(signal,l):
    x1=np.linspace(0,1,len(signal))
    
    piecewiseModel=Model(piecewise_prob)
    piecewiseModel.set_param_hint('l',     value=1,vary=False)
    piecewiseModel.set_param_hint('C',   vary=True,  value=.1,min=0,max=1)
    piecewiseModel.make_params()
    
    res=piecewiseModel.fit(signal,x=x1,weights=np.sin(1.5*x1)+1.5)
    return res


    

In [61]:

    

def save_config(lconf,conf):
    Config = configparser.ConfigParser()
    Config.optionxform = str
    
    cfgfile = open(join(lconf['node_path'],'project.conf'),'w')

    Config.add_section('Project')
    Config.set('Project','ProjectID',conf['project'])
    Config.set('Project','Cluster',conf['cluster'])
    Config.set('Project','JobName',conf['job_name'])
    Config.set('Project','JobNodes',conf['job_nodes'])
    Config.set('Project','CpuCores',conf['cpu_cores'])
    Config.set('Project','EstimatedTime',conf['estimated_time'])
    
    Config.add_section('Directories')
    Config.set('Directories','Node',conf['node_path'])
    Config.set('Directories','References',conf['ref_path'])
    Config.set('Directories','Data',conf['data_path'])
    Config.set('Directories','Output',conf['output_path'])
    Config.set('Directories','DoriC',conf['doric_path'])
    
    Config.add_section('Other')
    Config.set('Other','Mapper',conf['mapper'])
    Config.set('Other','RefName',conf['ref_name'])
    Config.set('Other','NrSamples',conf['nr_samples'])
    Config.set('Other','SamplesPerNode',conf['samples_per_node'])
    Config.set('Other','Email',conf['email'])
    Config.set('Other','DataPrefix',conf['data_prefix'])
    Config.set('Other','StartInd',conf['start_ind'])
    Config.set('Other','JobRange',conf['job_range'])
    Config.set('Other','FtpURL',conf['ftp_url'])
    Config.set('Other','DataURL',conf['data_url'])
    
    Config.write(cfgfile)
    cfgfile.close()


    

In [62]:

    

def mutation_run():
    rec.seq=mutateSeq(mutable_seq,.01)
    
    MutableSeq(str(rec.seq), IUPAC.unambiguous_dna)
    
    rec.id="mut_"+rec.id
    rec.description="mut_"+rec.description
    rec.name="mut_"+rec.name
    with open(sys.argv[2], "w") as output_handle:
        SeqIO.write(rec, output_handle, "fasta")


    

In [75]:

    

def compare_fit(c):
    err_hfit=[]
    err_pfit=[]
    res_fit=[]

    for i,samp in enumerate(c.samples):
        res_fit_tmp=[]
        err_hfit_tmp=[]
        err_pfit_tmp=[]
        
        for acc in c.pop.keys():
            try:
                depth_file=join(c.conf['output_path'],c.name+str(i),'npy',acc+'.depth.npy')
                best_file=join(c.conf['output_path'],c.name+str(i),'npy',acc+'.depth.best.npy')

                signal=2**(np.load(depth_file))
                signal=signal/(signal.sum()/len(signal))#*self.pop[acc].l)
                
                # extra deblurring
                signal=smooth_signal(signal,10**3)
                signal=sharpen_filter(signal)
                
                from_ptr=np.load(best_file)

                res=fit_signal(signal,c.pop[acc].l)
                res_fit_tmp.append(res.best_values['C'])

                err_hfit_tmp.append((c.pop[acc].C-res.best_values['C'])/c.pop[acc].C)

                err_pfit_tmp.append((c.pop[acc].C-from_ptr[2]+from_ptr[3])/c.pop[acc].C)
                #print "Simulated value: "+str(c.pop[acc].C)
                #print "Error from this fit: "+str(err_hfit_tmp[-1])+ ' value: ' + str(res.best_values['C'])
                #print "Error from initial PTR fit "+str(err_pfit_tmp[-1])+' value: ' + str(from_ptr[2]-from_ptr[3])
            except Exception as Ex:
                print Ex
                pass
        
        res_fit.append(res_fit_tmp)
        err_hfit.append(err_hfit_tmp)
        err_pfit.append(err_pfit_tmp)
    return [res_fit,np.abs(err_hfit),np.abs(err_pfit)]


    

In [64]:

    

def compare_fit_ptrc(c):
    ptrc=[]
    ptre=[]
    for i,samp in enumerate(c.samples):
        ptrc_tmp=[]
        ptre_tmp=[]
        cmd=['/bin/bash -c "']+['cd /mnt/vol ; /home/Programs/PTR-Pipeline/extra/ptrc_commands.sh']+[c.name+str(i)+'_1.fastq']+[c.name+str(i)+'_2.fastq"']
        create_mount_run('centos/hedani',td," ".join(cmd),envs)
        ptr_file=join(c.conf['output_path'],'ptrc_out')
        ptrc_frame=pd.read_pickle(ptr_file)
        ind_val=ptrc_frame.index.values
        for i in ind_val:#range(ptrc_frame.shape[0]):
            ptrc_tmp.append(ptrc_frame.loc[i]['sm'])#.ix[i]['sm'])
            ptre_tmp.append(np.abs(np.log2(ptrc_frame.loc[i]['sm'])-c.pop[i[:-1]].C)/c.pop[i[:-1]].C)
            
        ptrc.append(ptrc_tmp)
        ptre.append(ptre_tmp)
       
    return [np.log2(ptrc),ptre]


    

In [65]:

    

def mv_filt(L,omega):
    return (1/float(L))*(1-np.exp(-omega*1j*L))/(1-np.exp(-omega*1j))


    

In [66]:

    

def smooth_signal(signal,length):
    xi=np.linspace(0,1,len(signal))
    xi2=np.linspace(0,1,length)
    iss=scipy.interpolate.UnivariateSpline(xi,signal,k=3)
    iss.set_smoothing_factor(.00001)
    return iss(xi2)


    

In [67]:

    

def sharpen_filter(signal):
    x = np.linspace(-0.00001*np.pi/float(1), 1*np.pi/float(1),10**3)
    return np.roll(ifft(fft(signal,10**3)/mv_filt(1.27,x))[0:len(signal)],np.int64(-125/2))


    

In [68]:

    

def triangle(length, amplitude, y_offset):
    section = length // 4
    for direction in (1, -1):
        for i in range(section):
            yield y_offset + i * (amplitude / section) * direction
        for i in range(section):
            yield y_offset + (amplitude - (i * (amplitude / section))) * direction


    

In [ ]:

    




    

In [ ]:

    

df = pd.DataFrame([[1, 2], [3, 4]], columns=list('AB'))
df2 = pd.DataFrame([[5, 6], [7, 8]], columns=list('AB'))
df.append(df2)


    

In [ ]:

    

real_c=np.tile([growth_param[sel[0][1]][1],growth_param[sel[0][0]][1]],(5,1))
ptrc_e=np.abs(np.log2(ptrc)-real_c)/real_c


    

In [31]:

    

create_mount_run('centos/hedani',td," ".join(cmd),envs)


    

    
Out[31]:





'2016-12-11 20:40:13,175-ERROR: The following ids have been listed in the id_species files but not sequence was found: ti|821|org|NC_009614.1\nTraceback (most recent call last):\n  File "/home/Programs/PTRC/PTRC.py", line 138, in <module>\n    create_new_db(args.metadata, args.fasta, args.out)\n  File "cptr.py", line 629, in cptr.create_new_db (cptr.c:23369)\nAssertionError\nTraceback (most recent call last):\n  File "/home/Programs/PTRC/PTRC.py", line 133, in <module>\n    args.outfol, args.qual_format, args.cov_thresh, dbpath = args.db_path_name)\n  File "cptr.py", line 568, in cptr.coverage_analysis (cptr.c:20255)\n  File "cptr.py", line 226, in cptr._runGEM (cptr.c:9054)\n  File "cptr.py", line 185, in cptr._map_pe (cptr.c:7778)\n  File "cptr.py", line 203, in cptr._shell_command (cptr.c:8349)\n  File "/root/miniconda2/lib/python2.7/subprocess.py", line 574, in check_output\n    raise CalledProcessError(retcode, cmd, output=output)\nsubprocess.CalledProcessError: Command \'gem-mapper -I comm0/sim_db.gem -1 comm0/Data/comm00_1.fastq -2 comm0/Data/comm00_1.fastq -o comm0/Data/sm -q offset-33 -d all -p --max-extendable-matches all --max-extensions-per-match 5 -v\' returned non-zero exit status 1\nTraceback (most recent call last):\n  File "/home/Programs/PTRC/PTRC.py", line 136, in <module>\n    dbpath = args.db_path_name, csv_output = args.csv_output)\n  File "cptr.py", line 582, in cptr.calculate_ptrs (cptr.c:21032)\n  File "cptr.py", line 473, in cptr._load (cptr.c:16762)\n  File "/root/miniconda2/lib/python2.7/site-packages/pandas/io/pickle.py", line 65, in read_pickle\n    return try_read(path)\n  File "/root/miniconda2/lib/python2.7/site-packages/pandas/io/pickle.py", line 61, in try_read\n    with open(path, \'rb\') as fh:\nIOError: [Errno 2] No such file or directory: \'comm0/sim_db.pk\'\n'

In [747]:

    

cmd


    

    
Out[747]:





['/bin/bash -c "',
 'cd /mnt/vol ; /home/Programs/PTR-Pipeline/extra/ptrc_commands.sh',
 'comm04_1.fastq',
 'comm04_1.fastq"']

In [29]:

    

cmd=['/bin/bash -c "']+['cd /mnt/vol ; /home/Programs/PTR-Pipeline/extra/ptrc_commands.sh']+[c.name+str(0)+'_1.fastq']+[c.name+str(0)+'_1.fastq"']


    

In [30]:

    

" ".join(cmd)


    

    
Out[30]:





'/bin/bash -c " cd /mnt/vol ; /home/Programs/PTR-Pipeline/extra/ptrc_commands.sh comm00_1.fastq comm00_1.fastq"'

In [35]:

    

ptrc=pd.read_pickle(ptr_file)


    

In [36]:

    

ptrc.index.values


    

    
Out[36]:





array(['NC_000913.3\n'], dtype=object)

In [34]:

    

ptr_file=join(c.conf['output_path'],'ptrc_out')


    

In [716]:

    

ptrc.loc['Escherichia coli\n']['sm']


    

    
Out[716]:





1.3020408163265316

In [662]:

    

ptrc=compare_fit_ptrc(c)


    

In [685]:

    




    

In [69]:

    

acc=['NC_000913.3','NC_007779.1','NC_002655.2','NC_009614.1','NC_017218.1']
growth_param=[[.2, .5, .3, 100],[.1, .7, .2, 200],[.2, .8, .0001,250],[0.0001, .8, .2, 250],[.1, .7, .2,150]]
sel=[[0]]
td='/Users/Shared/CommTest'


    

In [70]:

    

cli = docker.Client(base_url='unix://var/run/docker.sock')
envs = docker.utils.parse_env_file(join(td,'env'))


    

In [71]:

    

totH=[]
totK=[]
for i in range(2):
    c=Community('comm0',[acc[s] for s in sel[0]],[growth_param[s] for s in sel[0]],td,'bowtie2',envs)

    c.build_index()

    tot_reads=np.linspace(5*10**3,10**4,2)
    for nr in tot_reads:
        c.sample(nr)

    c.write_reads()

    totK.append(compare_fit_ptrc(c))
    #[k_c,k_err]=

    c.run_pipeline()

    #[res,hfit,pfit]=
    
    
    totH.append(compare_fit(c))


    

In [85]:

    

totH[0]


    

    
Out[85]:





[[[[<lmfit.model.ModelResult at 0x1160a4ad0>],
   [<lmfit.model.ModelResult at 0x112b1d690>]],
  array([[ 0.05322417],
         [ 0.22688601]]),
  array([[ 0.21196545],
         [ 0.34804447]])],
 [[[<lmfit.model.ModelResult at 0x1160a4110>],
   [<lmfit.model.ModelResult at 0x1160afbd0>]],
  array([[ 0.03857808],
         [ 0.17221465]]),
  array([[ 0.17485312],
         [ 0.32419588]])]]

In [88]:

    

tips = sns.load_dataset("tips")


    

In [726]:

    

[k_c,k_err]=compare_fit_ptrc(c)


    

    




---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
<ipython-input-726-5d15a506d581> in <module>()
----> 1 [k_c,k_err]=compare_fit_ptrc(c)

<ipython-input-724-5b0971b9bed4> in compare_fit_ptrc(c)
      6         ptre_tmp=[]
      7         cmd=['/bin/bash -c "']+['cd /mnt/vol ; /home/Programs/PTR-Pipeline/extra/ptrc_commands.sh']+[c.name+str(i)+'_1.fastq']+[c.name+str(i)+'_2.fastq"']
----> 8         create_mount_run('centos/hedani',td," ".join(cmd),envs)
      9         ptr_file=join(c.conf['output_path'],'ptrc_out')
     10         ptrc_frame=pd.read_pickle(ptr_file)

<ipython-input-521-80ddf06e139f> in create_mount_run(image, mount_dir, cmd, envs)
     12     ctr=container.get('Id')
     13     cli.start(ctr)
---> 14     cli.wait(ctr,60*60*24*10)
     15     return cli.logs(ctr)

/Library/Python/2.7/site-packages/docker/utils/decorators.pyc in wrapped(self, resource_id, *args, **kwargs)
     19                 'image or container param is undefined'
     20             )
---> 21         return f(self, resource_id, *args, **kwargs)
     22     return wrapped
     23 

/Library/Python/2.7/site-packages/docker/api/container.pyc in wait(self, container, timeout)
    453     def wait(self, container, timeout=None):
    454         url = self._url("/containers/{0}/wait", container)
--> 455         res = self._post(url, timeout=timeout)
    456         self._raise_for_status(res)
    457         json_ = res.json()

/Library/Python/2.7/site-packages/docker/utils/decorators.pyc in inner(self, *args, **kwargs)
     45             else:
     46                 kwargs['headers'].update(self._auth_configs['HttpHeaders'])
---> 47         return f(self, *args, **kwargs)
     48     return inner

/Library/Python/2.7/site-packages/docker/client.pyc in _post(self, url, **kwargs)
    133     @update_headers
    134     def _post(self, url, **kwargs):
--> 135         return self.post(url, **self._set_request_timeout(kwargs))
    136 
    137     @update_headers

/Library/Python/2.7/site-packages/requests/sessions.pyc in post(self, url, data, json, **kwargs)
    533         """
    534 
--> 535         return self.request('POST', url, data=data, json=json, **kwargs)
    536 
    537     def put(self, url, data=None, **kwargs):

/Library/Python/2.7/site-packages/requests/sessions.pyc in request(self, method, url, params, data, headers, cookies, files, auth, timeout, allow_redirects, proxies, hooks, stream, verify, cert, json)
    486         }
    487         send_kwargs.update(settings)
--> 488         resp = self.send(prep, **send_kwargs)
    489 
    490         return resp

/Library/Python/2.7/site-packages/requests/sessions.pyc in send(self, request, **kwargs)
    607 
    608         # Send the request
--> 609         r = adapter.send(request, **kwargs)
    610 
    611         # Total elapsed time of the request (approximately)

/Library/Python/2.7/site-packages/requests/adapters.pyc in send(self, request, stream, timeout, verify, cert, proxies)
    421                     decode_content=False,
    422                     retries=self.max_retries,
--> 423                     timeout=timeout
    424                 )
    425 

/Library/Python/2.7/site-packages/requests/packages/urllib3/connectionpool.pyc in urlopen(self, method, url, body, headers, retries, redirect, assert_same_host, timeout, pool_timeout, release_conn, chunked, **response_kw)
    592                                                   timeout=timeout_obj,
    593                                                   body=body, headers=headers,
--> 594                                                   chunked=chunked)
    595 
    596             # If we're going to release the connection in ``finally:``, then

/Library/Python/2.7/site-packages/requests/packages/urllib3/connectionpool.pyc in _make_request(self, conn, method, url, timeout, chunked, **httplib_request_kw)
    382         try:
    383             try:  # Python 2.7, use buffering of HTTP responses
--> 384                 httplib_response = conn.getresponse(buffering=True)
    385             except TypeError:  # Python 2.6 and older, Python 3
    386                 try:

/usr/local/Cellar/python/2.7.11/Frameworks/Python.framework/Versions/2.7/lib/python2.7/httplib.pyc in getresponse(self, buffering)
   1134 
   1135         try:
-> 1136             response.begin()
   1137             assert response.will_close != _UNKNOWN
   1138             self.__state = _CS_IDLE

/usr/local/Cellar/python/2.7.11/Frameworks/Python.framework/Versions/2.7/lib/python2.7/httplib.pyc in begin(self)
    451         # read until we get a non-100 response
    452         while True:
--> 453             version, status, reason = self._read_status()
    454             if status != CONTINUE:
    455                 break

/usr/local/Cellar/python/2.7.11/Frameworks/Python.framework/Versions/2.7/lib/python2.7/httplib.pyc in _read_status(self)
    407     def _read_status(self):
    408         # Initialize with Simple-Response defaults
--> 409         line = self.fp.readline(_MAXLINE + 1)
    410         if len(line) > _MAXLINE:
    411             raise LineTooLong("header line")

/usr/local/Cellar/python/2.7.11/Frameworks/Python.framework/Versions/2.7/lib/python2.7/socket.pyc in readline(self, size)
    478             while True:
    479                 try:
--> 480                     data = self._sock.recv(self._rbufsize)
    481                 except error, e:
    482                     if e.args[0] == EINTR:

KeyboardInterrupt: 

In [708]:

    

plt.plot(tot_reads,np.mean(hfit,1),tot_reads,np.mean(ptrc_e,1))#,tot_reads,[r[0].best_values['C'],r[1].best_values['C'] for r in res])
legend = plt.legend(['H','K'])
#ax = plt.axes()
#ax.set_title('Mean absolute error')
#g.fig.suptitle('THIS IS A TITLE, YOU BET')


    

In [ ]:

    

sns_plot.savefig("output.pdf")


    

In [280]:

    

args0=['/home/Programs/PTR-Pipeline/ptr_pipeline.py','-c',join('/mnt/vol','comm0','project.conf'),'-e','hedani@chalmers.se']          
create_mount_run('centos/hedani',td," ".join(['/bin/bash -c "']+args0+['make']+[';']+args0+['run"']),envs)


    

    
Out[280]:





"Running alignment stage\ncomm00:\nAligning comm00:\n99999 reads; of these:\n  99999 (100.00%) were paired; of these:\n    99690 (99.69%) aligned concordantly 0 times\n    147 (0.15%) aligned concordantly exactly 1 time\n    162 (0.16%) aligned concordantly >1 times\n    ----\n    99690 pairs aligned concordantly 0 times; of these:\n      47312 (47.46%) aligned discordantly 1 time\n    ----\n    52378 pairs aligned 0 times concordantly or discordantly; of these:\n      104756 mates make up the pairs; of these:\n        57 (0.05%) aligned 0 times\n        4319 (4.12%) aligned exactly 1 time\n        100380 (95.82%) aligned >1 times\n99.97% overall alignment rate\nMoving into subjobs\ncomm00: Running pathoscope stage\nRunning without mySQLdb library\ncomm00: Running coverage stage\ncomm00: Sorting\ncomm00: Calculating coverage\ncomm00: Performing piecewise fits\nNC_000913.3\nNC_000913.3.depth\nNC_000913.3.depth: Bin coverage too low, exiting.\nNC_009614.1\nNC_009614.1.depth\nNC_009614.1.depth: Coverage OK (0.813975626304 x).\n\ncntr: 6\nNC_002655.2\nNC_002655.2.depth\n/home/Programs/PTR-Pipeline/bin/piecewiseFit.py:65: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future\n  out=np.zeros(iEnd);\n/home/Programs/PTR-Pipeline/bin/piecewiseFit.py:72: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future\n  out=np.zeros(np.floor(len(x)/slide));\n/home/Programs/PTR-Pipeline/bin/piecewiseFit.py:65: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future\n  out=np.zeros(iEnd);\n/home/Programs/PTR-Pipeline/bin/piecewiseFit.py:72: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future\n  out=np.zeros(np.floor(len(x)/slide));\nmv: cannot stat 'depth/*.log': No such file or directory\nNC_002655.2.depth: Coverage OK (0.810297289744 x).\n\ncntr: 8\ncomm00: Performing secondary fits\n*\nI/O error *.npy\n"

In [507]:

    

res=fit_signal(st,1)
err_hfit=((.5-res.best_values['C'])/.5)
err_hfit


    

    
Out[507]:





0.8969818062406385

In [464]:

    

res.best_values['C']


    

    
Out[464]:





0.50223664098758436

In [496]:

    

#depth_file=join(c.conf['output_path'],c.name+'0','npy',c.pop.keys()[0]+'.depth.npy')
depth_file='/Users/Shared/CommTest/comm0/Out/comm00/npy/NC_000913.3.depth.npy'


    

In [508]:

    

signal=2**np.load(depth_file)
signal=signal/(signal.sum()/len(signal))


    

In [509]:

    

plt.plot(signal)


    

    
Out[509]:





[<matplotlib.lines.Line2D at 0x137fd01d0>]






    

In [429]:

    

x = np.linspace(-2*np.pi, 2*np.pi, len(signal))


    

In [339]:

    

xx=np.linspace(-np.pi, np.pi, 1000)
plt.plot(np.abs(mv_filt(1.2,xx)))


    

    
Out[339]:





[<matplotlib.lines.Line2D at 0x12cf47f10>]






    

In [401]:

    

y2=[]
x2=np.linspace(0,.5,100)
y2=np.append(pxn2(x2,.5),pxn2(x2,.5)[::-1])
x2=np.linspace(0,1,200)

y=list(triangle(1000, 0.1,10))
x=range(len(y))
x_p = np.linspace(-2*np.pi, 2*np.pi, len(y2))

f=mv_filt(20,x_p)
fft_y=fft(y2)*f
smooth_y=np.roll(ifft(fft_y),0)
smooth_y=smooth_y/np.sum(smooth_y)

x3=np.linspace(0,1,len(y_long))
y_skum=y_long/(np.sum(y_long+np.mean(y_long)))
x_p2 = np.linspace(-2*np.pi, 2*np.pi, len(y_skum))

plt.plot(x2,y2/np.sum(y2),x2,smooth_y,x2,ifft(fft(smooth_y)/f))


    

    
Out[401]:





[<matplotlib.lines.Line2D at 0x1325a21d0>,
 <matplotlib.lines.Line2D at 0x1325a2490>,
 <matplotlib.lines.Line2D at 0x1325a2b50>]






    

In [505]:

    

xi=np.linspace(0,1,len(signal))
xi2=np.linspace(0,1,10**3)
iss=scipy.interpolate.UnivariateSpline(xi,signal,k=3)
iss.set_smoothing_factor(.00001)
y_long=iss(xi2)
plt.plot(y_long)


    

    
Out[505]:





[<matplotlib.lines.Line2D at 0x137ce6410>]






    

In [208]:

    

plt.plot(np.abs(ifft(fft(signal,10**4))[0:len(signal)]))


    

    
Out[208]:





[<matplotlib.lines.Line2D at 0x121a0eed0>]






    

In [512]:

    

signal=y_long
x = np.linspace(-0.00001*np.pi/float(1), 1*np.pi/float(1),10**3)
st=np.roll(ifft(fft(signal,10**3)/mv_filt(1.27,x))[0:len(signal)],np.int64(-125/2))
x2=np.linspace(0,1,len(st))
x1=np.linspace(0,1,len(signal))
plt.plot(x2,st/np.sum(st),x1,signal/np.sum(signal))


    

    
Out[512]:





[<matplotlib.lines.Line2D at 0x1383f14d0>,
 <matplotlib.lines.Line2D at 0x1383f1790>]






    

In [514]:

    

plt.plot(x2,st/np.sum(st),x1,signal/np.sum(signal))


    

    
Out[514]:





[<matplotlib.lines.Line2D at 0x1384a5c10>,
 <matplotlib.lines.Line2D at 0x1384a5ed0>]






    

In [193]:

    

plt.plot(np.abs(mv_filt(1*10**1,x)))


    

    
Out[193]:





[<matplotlib.lines.Line2D at 0x12058b790>]






    

In [411]:

    

np.sqrt(-1)


    

    
Out[411]:





nan

In [ ]:

    

1+np.exp()


    

In [398]:

    

np.exp(1)


    

    
Out[398]:





2.7182818284590451

In [272]:

    

create_mount_run('centos/hedani',td," ".join(['/bin/bash -c "']+['which python"']),envs)


    

    
Out[272]:





'/usr/bin/python\n'

In [52]:

    

tot_reads=[10**5]#np.linspace(5*10**3,10**4,15)
c=Community('comm0',[acc[s] for s in sel[0]],[growth_param[s] for s in sel[0]],td,'bowtie2')
c.build_index()
for nr in tot_reads:
    c.sample(nr)
    c.write_reads()

c.run_pipeline()
[res,hfit,pfit]=c.compare_fit()


    

    




---------------------------------------------------------------------------
NotFound                                  Traceback (most recent call last)
<ipython-input-52-071c0a386248> in <module>()
      1 tot_reads=[10**5]#np.linspace(5*10**3,10**4,15)
----> 2 c=Community('comm0',[acc[s] for s in sel[0]],[growth_param[s] for s in sel[0]],td,'bowtie2')
      3 c.build_index()
      4 for nr in tot_reads:
      5     c.sample(nr)

<ipython-input-31-7a930caef3da> in __init__(self, name, acc, growth_param, td, mapper)
      6 
      7         self.create_dirs(True)
----> 8         self.fetch_ref(acc)
      9 
     10         self.pop = self.init_population(acc,growth_param)

<ipython-input-31-7a930caef3da> in fetch_ref(self, acc)
    100         f.close()
    101 
--> 102         create_mount_run('centos/hedani','/Users/Shared/d_py',['-e','hedani@chalmers.se','fetch-references','-s',acc_path])
    103         #args=ptr_pipeline.parse_args(['-e','hedani@chalmers.se','fetch-references','-s',acc_path])
    104         #ptr_pipeline.main(args,self.conf)

<ipython-input-28-257136503bbf> in create_mount_run(image, mount_dir, cmd)
      9         })
     10     )
---> 11     cli.start(container=container.get('Id'))
     12     return cli.logs(container=container.get('Id'))

/Library/Python/2.7/site-packages/docker/utils/decorators.pyc in wrapped(self, resource_id, *args, **kwargs)
     19                 'image or container param is undefined'
     20             )
---> 21         return f(self, resource_id, *args, **kwargs)
     22     return wrapped
     23 

/Library/Python/2.7/site-packages/docker/api/container.pyc in start(self, container, binds, port_bindings, lxc_conf, publish_all_ports, links, privileged, dns, dns_search, volumes_from, network_mode, restart_policy, cap_add, cap_drop, devices, extra_hosts, read_only, pid_mode, ipc_mode, security_opt, ulimits)
    381         url = self._url("/containers/{0}/start", container)
    382         res = self._post_json(url, data=start_config)
--> 383         self._raise_for_status(res)
    384 
    385     @utils.minimum_version('1.17')

/Library/Python/2.7/site-packages/docker/client.pyc in _raise_for_status(self, response, explanation)
    171         except requests.exceptions.HTTPError as e:
    172             if e.response.status_code == 404:
--> 173                 raise errors.NotFound(e, response, explanation=explanation)
    174             raise errors.APIError(e, response, explanation=explanation)
    175 

NotFound: 404 Client Error: Not Found ("{"message":"oci runtime error: exec: \"-e\": executable file not found in $PATH"}")

In [ ]:

    




    

In [247]:

    

community=[]
nrs=np.linspace(10**5,10**5,1)
res=[]
for nr in nrs:
    community.append(run_community_analysis('fs',td,'bowtie2',acc,growth_param,sel,nr,0))
    #res.append(community[i][0].compare_fit())


    

    




ERROR:root:Internal Python error in the inspect module.
Below is the traceback from this internal error.







    




Traceback (most recent call last):
  File "/Library/Python/2.7/site-packages/IPython/core/ultratb.py", line 1132, in get_records
    return _fixed_getinnerframes(etb, number_of_lines_of_context, tb_offset)
  File "/Library/Python/2.7/site-packages/IPython/core/ultratb.py", line 313, in wrapped
    return f(*args, **kwargs)
  File "/Library/Python/2.7/site-packages/IPython/core/ultratb.py", line 358, in _fixed_getinnerframes
    records = fix_frame_records_filenames(inspect.getinnerframes(etb, context))
  File "/usr/local/Cellar/python/2.7.11/Frameworks/Python.framework/Versions/2.7/lib/python2.7/inspect.py", line 1049, in getinnerframes
    framelist.append((tb.tb_frame,) + getframeinfo(tb, context))
  File "/usr/local/Cellar/python/2.7.11/Frameworks/Python.framework/Versions/2.7/lib/python2.7/inspect.py", line 1009, in getframeinfo
    filename = getsourcefile(frame) or getfile(frame)
  File "/usr/local/Cellar/python/2.7.11/Frameworks/Python.framework/Versions/2.7/lib/python2.7/inspect.py", line 454, in getsourcefile
    if hasattr(getmodule(object, filename), '__loader__'):
  File "/usr/local/Cellar/python/2.7.11/Frameworks/Python.framework/Versions/2.7/lib/python2.7/inspect.py", line 483, in getmodule
    file = getabsfile(object, _filename)
  File "/usr/local/Cellar/python/2.7.11/Frameworks/Python.framework/Versions/2.7/lib/python2.7/inspect.py", line 467, in getabsfile
    return os.path.normcase(os.path.abspath(_filename))
  File "/usr/local/Cellar/python/2.7.11/Frameworks/Python.framework/Versions/2.7/lib/python2.7/posixpath.py", line 364, in abspath
    cwd = os.getcwd()
OSError: [Errno 2] No such file or directory






    




---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
/Library/Python/2.7/site-packages/IPython/core/interactiveshell.pyc in run_code(self, code_obj, result)
   2896             if result is not None:
   2897                 result.error_in_exec = sys.exc_info()[1]
-> 2898             self.showtraceback()
   2899         else:
   2900             outflag = 0

/Library/Python/2.7/site-packages/IPython/core/interactiveshell.pyc in showtraceback(self, exc_tuple, filename, tb_offset, exception_only)
   1822                     except Exception:
   1823                         stb = self.InteractiveTB.structured_traceback(etype,
-> 1824                                             value, tb, tb_offset=tb_offset)
   1825 
   1826                     self._showtraceback(etype, value, stb)

/Library/Python/2.7/site-packages/IPython/core/ultratb.pyc in structured_traceback(self, etype, value, tb, tb_offset, number_of_lines_of_context)
   1404         self.tb = tb
   1405         return FormattedTB.structured_traceback(
-> 1406             self, etype, value, tb, tb_offset, number_of_lines_of_context)
   1407 
   1408 

/Library/Python/2.7/site-packages/IPython/core/ultratb.pyc in structured_traceback(self, etype, value, tb, tb_offset, number_of_lines_of_context)
   1312             # Verbose modes need a full traceback
   1313             return VerboseTB.structured_traceback(
-> 1314                 self, etype, value, tb, tb_offset, number_of_lines_of_context
   1315             )
   1316         else:

/Library/Python/2.7/site-packages/IPython/core/ultratb.pyc in structured_traceback(self, etype, evalue, etb, tb_offset, number_of_lines_of_context)
   1196                 structured_traceback_parts += formatted_exception
   1197         else:
-> 1198             structured_traceback_parts += formatted_exception[0]
   1199 
   1200         return structured_traceback_parts

IndexError: string index out of range

In [73]:

    

c=Community('comm0',[acc[s] for s in sel[0]],[growth_param[s] for s in sel[0]],td,'bowtie2')


    

    




---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
<ipython-input-73-fde9aa008fdc> in <module>()
----> 1 c=Community('comm0',[acc[s] for s in sel[0]],[growth_param[s] for s in sel[0]],td,'bowtie2')

<ipython-input-55-d175fa22a917> in __init__(self, name, acc, growth_param, td, mapper)
     12         self.fetch_ref(acc)
     13 
---> 14         self.pop = self.init_population(acc,growth_param)
     15         self.distribution = self.community_distribution()
     16         self.samples = []

<ipython-input-55-d175fa22a917> in init_population(self, acc, growth_param)
     29         for i,a in enumerate(acc):
     30             pop[a]=population(B = growth_param[i][0],C = growth_param[i][1],
---> 31                              D = growth_param[i][2], l = len(records[a]),
     32                              seq = records[a], cells=growth_param[i][3])
     33         return pop

KeyError: 'NC_000913.3'

In [47]:

    

tot_reads=[10**5]#np.linspace(5*10**3,10**4,15)
c=Community('comm0',[acc[s] for s in sel[0]],[growth_param[s] for s in sel[0]],td,'bowtie2')
c.build_index()


    

    




---------------------------------------------------------------------------
IOError                                   Traceback (most recent call last)
<ipython-input-47-5b05c90de0d8> in <module>()
      1 tot_reads=[10**5]#np.linspace(5*10**3,10**4,15)
----> 2 c=Community('comm0',[acc[s] for s in sel[0]],[growth_param[s] for s in sel[0]],td,'bowtie2')
      3 c.build_index()

<ipython-input-29-c0505f730c20> in __init__(self, name, acc, growth_param, td, mapper)
      7         self.args0=['/home/Programs/PTR-Pipeline/ptr_pipeline.py','-c',join(self.d_conf['node_path'],'project.conf'),'-e','hedani@chalmers.se']
      8 
----> 9         save_config(self.d_conf,self.conf)
     10 
     11         self.create_dirs(True)

<ipython-input-43-68021bfe077a> in save_config(lconf, conf)
      3     Config.optionxform = str
      4 
----> 5     cfgfile = open(join(lconf['node_path'],'project.conf'),'w')
      6 
      7     Config.add_section('Project')

IOError: [Errno 2] No such file or directory: '/mnt/vol/comm0/project.conf'

In [2132]:

    

n=130;
D1_snok=np.random.randint(30,40);C1_snok=np.random.randint(40,70);tau1_snok=np.random.randint(130,250,n)
#tau1_snok[-1]=50;
D2_snok=10.;C2_snok=35.;tau2_snok=np.random.randint(60,200,n)
ab1=[np.random.rand() for x in range(n)];ab2=np.random.randint(100,200,n);
org1=[[np.max([1-(C1_snok+D1_snok)/float(tau1_snok[i]),0]),(C1_snok)/float(tau1_snok[i]), (D1_snok)/float(tau1_snok[i]),ab1[i]] for i in range(tau1_snok.size)]
org2=[[1-(C2_snok+D2_snok)/float(tau2_snok[i]),(C2_snok)/float(tau2_snok[i]), (D1_snok)/float(tau2_snok[i]),ab2[i]/float(ab1[i])] for i in range(tau2_snok.size)]

#same_org[0][0]=0


    

In [2133]:

    

print D1_snok,C1_snok


    

    




34 44

In [2134]:

    

same_org=np.concatenate((org1,org2))


    

In [2135]:

    

Ri1


    

    
Out[2135]:





array([ 0.12517895,  0.28003121,  0.11440445,  0.15001599,  0.07173145,
        0.32649918,  0.39299506,  0.01202917,  0.34304507,  0.33282179,
        0.04329798,  0.26181841,  0.33872681,  0.14846477,  0.24897509,
        0.061663  ,  0.2401119 ,  0.40135043,  0.2886033 ,  0.39829411,
        0.1850574 ,  0.34708688,  0.20622304,  0.16787013,  0.16767945,
        0.37298813,  0.33103408,  0.25134168,  0.35186729,  0.16303172,
        0.30711971,  0.26330621,  0.12609793,  0.20625785,  0.3764792 ,
        0.22850401,  0.0450059 ,  0.0809682 ,  0.18655135,  0.00784618,
        0.29397322,  0.08577948,  0.29319743,  0.10284626,  0.30304625,
        0.28011791,  0.1960176 ,  0.16534839,  0.02916867,  0.0991715 ,
        0.0578688 ,  0.27499659,  0.027967  ,  0.13073409,  0.15571648,
        0.01069903,  0.1112608 ,  0.18262335,  0.20584184,  0.3887354 ,
        0.38976375,  0.38220299,  0.36039617,  0.10669672,  0.22368446,
        0.23371822,  0.03268177,  0.25700886,  0.07782686,  0.37210099,
        0.11653302,  0.06587968,  0.14566878,  0.0068028 ,  0.22986611,
        0.00415846,  0.21185567,  0.23742862,  0.36762303,  0.07634184,
        0.11859662,  0.3574003 ,  0.13858283,  0.08572636,  0.01833882,
        0.19050518,  0.34799022,  0.24642473,  0.10322462,  0.31964564,
        0.07285449,  0.36825552,  0.09071813,  0.18713142,  0.23090128,
        0.19523249,  0.24599339,  0.2281989 ,  0.0326502 ,  0.04275451,
        0.29883877,  0.34733176,  0.25028025,  0.24900968,  0.01261928,
        0.3568365 ,  0.01669861,  0.29840329,  0.10386471,  0.27002372,
        0.04221017,  0.23762403,  0.08465263,  0.40777278,  0.09767022,
        0.12452879,  0.31985767,  0.09076833,  0.15271323,  0.30922279,
        0.2969826 ,  0.29776965,  0.21090022,  0.3868555 ,  0.17569115,
        0.27928653,  0.43143433,  0.0779982 ,  0.01669648,  0.10084281])

In [2136]:

    

Ri1=np.array([Gekv2(p[1],p[2])*1*p[3] for p in org1])
rs=np.random.randint(2,3,n);
for i in range(n):
    Ri1[i]=Ri1[i]/rs[i];
    org1[i][3]=org1[i][3]/rs[i];


Ki1=np.array([(2**(p[1])-1)/(float(p[1])*np.log(2)) for p in org1])
Ri2=np.array([Gekv2(p[1],p[2])*1*p[3] for p in org2])
Ki2=np.array([(2**(p[1])-1)/(float(p[1])*np.log(2)) for p in org2])


    

In [2137]:

    

org1


    

    
Out[2137]:





[[0.64864864864864868,
  0.1981981981981982,
  0.15315315315315314,
  0.46079501691747554],
 [0.59375, 0.22916666666666666, 0.17708333333333334, 0.20116365324132957],
 [0.67088607594936711,
  0.18565400843881857,
  0.14345991561181434,
  0.49654154079404594],
 [0.57837837837837835,
  0.23783783783783785,
  0.1837837837837838,
  0.24764121123035404],
 [0.60999999999999999, 0.22, 0.17, 0.3921142465319019],
 [0.6637931034482758,
  0.1896551724137931,
  0.14655172413793102,
  0.20158666914294182],
 [0.6637931034482758,
  0.1896551724137931,
  0.14655172413793102,
  0.22483740954130388],
 [0.65638766519823788,
  0.19383259911894274,
  0.14977973568281938,
  0.22464293403338875],
 [0.50943396226415094,
  0.27672955974842767,
  0.2138364779874214,
  0.44063274576640726],
 [0.60204081632653061,
  0.22448979591836735,
  0.17346938775510204,
  0.43008967519542796],
 [0.60606060606060608,
  0.2222222222222222,
  0.1717171717171717,
  0.17211265559073496],
 [0.52439024390243905,
  0.2682926829268293,
  0.2073170731707317,
  0.2866343753302563],
 [0.58288770053475936,
  0.23529411764705882,
  0.18181818181818182,
  0.17988289772514116],
 [0.45454545454545459,
  0.3076923076923077,
  0.23776223776223776,
  0.0918123170900288],
 [0.45454545454545459,
  0.3076923076923077,
  0.23776223776223776,
  0.4202347320171348],
 [0.58730158730158732,
  0.2328042328042328,
  0.17989417989417988,
  0.4089650647473905],
 [0.59999999999999998,
  0.22564102564102564,
  0.17435897435897435,
  0.4737860668884134],
 [0.68163265306122445,
  0.17959183673469387,
  0.13877551020408163,
  0.18682768873395322],
 [0.53846153846153844,
  0.2603550295857988,
  0.20118343195266272,
  0.46915131141688116],
 [0.58947368421052637,
  0.23157894736842105,
  0.17894736842105263,
  0.14030413639900796],
 [0.53012048192771077,
  0.26506024096385544,
  0.20481927710843373,
  0.01801941652267519],
 [0.64220183486238525,
  0.2018348623853211,
  0.1559633027522936,
  0.46220745532265606],
 [0.64055299539170507,
  0.20276497695852536,
  0.15668202764976957,
  0.11793696321699298],
 [0.40000000000000002,
  0.3384615384615385,
  0.26153846153846155,
  0.016939443057983616],
 [0.46575342465753422,
  0.3013698630136986,
  0.2328767123287671,
  0.02568227640343973],
 [0.43065693430656937,
  0.32116788321167883,
  0.24817518248175183,
  0.22051022260961034],
 [0.60406091370558368,
  0.2233502538071066,
  0.17258883248730963,
  0.012433187891947306],
 [0.6080402010050252,
  0.22110552763819097,
  0.1708542713567839,
  0.2914623529141575],
 [0.46938775510204078,
  0.29931972789115646,
  0.23129251700680273,
  0.0003985532846665918],
 [0.61194029850746268,
  0.21890547263681592,
  0.1691542288557214,
  0.4267106671672934],
 [0.67901234567901236,
  0.18106995884773663,
  0.13991769547325103,
  0.4075861771070389],
 [0.67901234567901236,
  0.18106995884773663,
  0.13991769547325103,
  0.17460659919835086],
 [0.45454545454545459,
  0.3076923076923077,
  0.23776223776223776,
  0.23454114179864116],
 [0.63888888888888884,
  0.2037037037037037,
  0.1574074074074074,
  0.07579110856008303],
 [0.5357142857142857,
  0.2619047619047619,
  0.20238095238095238,
  0.0652801694669164],
 [0.66808510638297869,
  0.18723404255319148,
  0.14468085106382977,
  0.27763616926454493],
 [0.68674698795180722,
  0.17670682730923695,
  0.13654618473895583,
  0.022855740658188184],
 [0.58288770053475936,
  0.23529411764705882,
  0.18181818181818182,
  0.49074281480303383],
 [0.65638766519823788,
  0.19383259911894274,
  0.14977973568281938,
  0.02739239734518567],
 [0.45454545454545459,
  0.3076923076923077,
  0.23776223776223776,
  0.4917295420246727],
 [0.66233766233766234,
  0.19047619047619047,
  0.1471861471861472,
  0.2497526462054266],
 [0.6080402010050252,
  0.22110552763819097,
  0.1708542713567839,
  0.4770605142372987],
 [0.43478260869565222,
  0.3188405797101449,
  0.2463768115942029,
  0.28254111584207114],
 [0.45833333333333337,
  0.3055555555555556,
  0.2361111111111111,
  0.010230785310108625],
 [0.61576354679802958,
  0.21674876847290642,
  0.16748768472906403,
  0.291575700978587],
 [0.62318840579710144,
  0.21256038647342995,
  0.1642512077294686,
  0.1643615088854778],
 [0.47651006711409394,
  0.2953020134228188,
  0.22818791946308725,
  0.285877263586],
 [0.4135338345864662,
  0.3308270676691729,
  0.2556390977443609,
  0.4079787945674634],
 [0.49677419354838714,
  0.2838709677419355,
  0.21935483870967742,
  0.21611641119296693],
 [0.59999999999999998,
  0.22564102564102564,
  0.17435897435897435,
  0.31412258182920005],
 [0.58288770053475936,
  0.23529411764705882,
  0.18181818181818182,
  0.45168765515216625],
 [0.68421052631578949,
  0.17813765182186234,
  0.13765182186234817,
  0.2223351304214118],
 [0.46575342465753422,
  0.3013698630136986,
  0.2328767123287671,
  0.06362377612636472],
 [0.68032786885245899,
  0.18032786885245902,
  0.13934426229508196,
  0.31736353876696866],
 [0.50632911392405067,
  0.27848101265822783,
  0.21518987341772153,
  0.011060227802662803],
 [0.63207547169811318,
  0.20754716981132076,
  0.16037735849056603,
  0.40062911535554124],
 [0.62318840579710144,
  0.21256038647342995,
  0.1642512077294686,
  0.044284671774929896],
 [0.62318840579710144,
  0.21256038647342995,
  0.1642512077294686,
  0.017877932941876062],
 [0.4388489208633094,
  0.31654676258992803,
  0.2446043165467626,
  0.27840805982487965],
 [0.46575342465753422,
  0.3013698630136986,
  0.2328767123287671,
  0.4303743176937124],
 [0.57608695652173914,
  0.2391304347826087,
  0.18478260869565216,
  0.4789520760689538],
 [0.62857142857142856,
  0.20952380952380953,
  0.1619047619047619,
  0.39675886455683007],
 [0.58288770053475936,
  0.23529411764705882,
  0.18181818181818182,
  0.17842565107350478],
 [0.45454545454545459,
  0.3076923076923077,
  0.23776223776223776,
  0.15424981483099498],
 [0.68292682926829262,
  0.17886178861788618,
  0.13821138211382114,
  0.4159795866799625],
 [0.47297297297297303,
  0.2972972972972973,
  0.22972972972972974,
  0.0400130749241025],
 [0.67088607594936711,
  0.18565400843881857,
  0.14345991561181434,
  0.16761038698281472],
 [0.55932203389830515,
  0.24858757062146894,
  0.192090395480226,
  0.06748110497628335],
 [0.50943396226415094,
  0.27672955974842767,
  0.2138364779874214,
  0.33515107188402804],
 [0.53012048192771077,
  0.26506024096385544,
  0.20481927710843373,
  0.37901848137513133],
 [0.64864864864864868,
  0.1981981981981982,
  0.15315315315315314,
  0.32363978989470665],
 [0.64055299539170507,
  0.20276497695852536,
  0.15668202764976957,
  0.32810332357826466],
 [0.52727272727272734,
  0.26666666666666666,
  0.20606060606060606,
  0.4715145827317281],
 [0.56666666666666665,
  0.24444444444444444,
  0.18888888888888888,
  0.004128507612861987],
 [0.61764705882352944,
  0.21568627450980393,
  0.16666666666666666,
  0.25827871302554967],
 [0.67226890756302526,
  0.18487394957983194,
  0.14285714285714285,
  0.4987618526306502],
 [0.6637931034482758,
  0.1896551724137931,
  0.14655172413793102,
  0.47396175599898605],
 [0.67901234567901236,
  0.18106995884773663,
  0.13991769547325103,
  0.22916006140325845],
 [0.66808510638297869,
  0.18723404255319148,
  0.14468085106382977,
  0.006042670814164486],
 [0.61764705882352944,
  0.21568627450980393,
  0.16666666666666666,
  0.31867207806094916],
 [0.63033175355450233,
  0.20853080568720378,
  0.16113744075829384,
  0.4533736478201163],
 [0.56906077348066297,
  0.2430939226519337,
  0.1878453038674033,
  0.2480520511184005],
 [0.68163265306122445,
  0.17959183673469387,
  0.13877551020408163,
  0.16669709252320397],
 [0.57608695652173914,
  0.2391304347826087,
  0.18478260869565216,
  0.46919015186765267],
 [0.63551401869158886,
  0.205607476635514,
  0.1588785046728972,
  0.380183644923833],
 [0.64864864864864868,
  0.1981981981981982,
  0.15315315315315314,
  0.2307091141351154],
 [0.58064516129032251,
  0.23655913978494625,
  0.1827956989247312,
  0.16310104203088738],
 [0.55681818181818188, 0.25, 0.19318181818181818, 0.29090402131285625],
 [0.55932203389830515,
  0.24858757062146894,
  0.192090395480226,
  0.3036843665003522],
 [0.66949152542372881,
  0.1864406779661017,
  0.1440677966101695,
  0.3260796666861614],
 [0.67901234567901236,
  0.18106995884773663,
  0.13991769547325103,
  0.25079441129544927],
 [0.46938775510204078,
  0.29931972789115646,
  0.23129251700680273,
  0.05843929831879591],
 [0.56906077348066297,
  0.2430939226519337,
  0.1878453038674033,
  0.3788527812958234],
 [0.51552795031055898,
  0.2732919254658385,
  0.2111801242236025,
  0.2422078913952499],
 [0.60999999999999999, 0.22, 0.17, 0.4672477014175954],
 [0.56666666666666665,
  0.24444444444444444,
  0.18888888888888888,
  0.42521132206742074],
 [0.50943396226415094,
  0.27672955974842767,
  0.2138364779874214,
  0.10407675396895005],
 [0.58947368421052637,
  0.23157894736842105,
  0.17894736842105263,
  0.47486698927718024],
 [0.68163265306122445,
  0.17959183673469387,
  0.13877551020408163,
  0.07495626881189427],
 [0.42647058823529416, 0.3235294117647059, 0.25, 0.48954133636693964],
 [0.54117647058823537, 0.25882352941176473, 0.2, 0.3871470533255098],
 [0.46575342465753422,
  0.3013698630136986,
  0.2328767123287671,
  0.2675470437722781],
 [0.65333333333333332,
  0.19555555555555557,
  0.1511111111111111,
  0.3255599590818351],
 [0.64864864864864868,
  0.1981981981981982,
  0.15315315315315314,
  0.45573537792144153],
 [0.63720930232558137,
  0.20465116279069767,
  0.15813953488372093,
  0.17392641682893523],
 [0.6637931034482758,
  0.1896551724137931,
  0.14655172413793102,
  0.2968118195696663],
 [0.65333333333333332,
  0.19555555555555557,
  0.1511111111111111,
  0.07989447253978554],
 [0.51552795031055898,
  0.2732919254658385,
  0.2111801242236025,
  0.49371289448667205],
 [0.47999999999999998,
  0.29333333333333333,
  0.22666666666666666,
  0.09733296432041799],
 [0.5185185185185186,
  0.2716049382716049,
  0.20987654320987653,
  0.26059767520839067],
 [0.46938775510204078,
  0.29931972789115646,
  0.23129251700680273,
  0.14626166795275647],
 [0.66233766233766234,
  0.19047619047619047,
  0.1471861471861472,
  0.1628807839145976],
 [0.60999999999999999, 0.22, 0.17, 0.16708772481168682],
 [0.65789473684210531,
  0.19298245614035087,
  0.14912280701754385,
  0.26328747680426245],
 [0.68674698795180722,
  0.17670682730923695,
  0.13654618473895583,
  0.1593009524224957],
 [0.56424581005586594,
  0.24581005586592178,
  0.18994413407821228,
  0.460562307244909],
 [0.61386138613861385,
  0.21782178217821782,
  0.16831683168316833,
  0.28121643054034673],
 [0.48684210526315785,
  0.2894736842105263,
  0.2236842105263158,
  0.1253352606794409],
 [0.45833333333333337,
  0.3055555555555556,
  0.2361111111111111,
  0.010656670378031252],
 [0.56666666666666665,
  0.24444444444444444,
  0.18888888888888888,
  0.44690907408942854],
 [0.53012048192771077,
  0.26506024096385544,
  0.20481927710843373,
  0.31841682872717114],
 [0.65333333333333332,
  0.19555555555555557,
  0.1511111111111111,
  0.4744375311911005],
 [0.68292682926829262,
  0.17886178861788618,
  0.13821138211382114,
  0.31763678516257166],
 [0.49350649350649356,
  0.2857142857142857,
  0.22077922077922077,
  0.3042811363098526],
 [0.61951219512195121,
  0.2146341463414634,
  0.16585365853658537,
  0.10414107380361959],
 [0.40000000000000002,
  0.3384615384615385,
  0.26153846153846155,
  0.4481668632525234],
 [0.55932203389830515,
  0.24858757062146894,
  0.192090395480226,
  0.27911572937343504],
 [0.60406091370558368,
  0.2233502538071066,
  0.17258883248730963,
  0.19864426249618172],
 [0.50943396226415094,
  0.27672955974842767,
  0.2138364779874214,
  0.19862820800678094],
 [0.52727272727272734,
  0.26666666666666666,
  0.20606060606060606,
  0.22704483674014786]]

In [1887]:

    

np.array([Gekv2(p[1],p[2])*1*p[3] for p in org1])


    

    
Out[1887]:





array([ 1.09875756,  0.54429242,  0.14761133])

In [1888]:

    

[2**(p[2])*((2**(p[1])-1)/(float(p[1])*np.log(2)))*1*p[3]/400. for p in same_org]


    

    
Out[1888]:





[0.002746893895910563,
 0.0013607310513527526,
 0.00036902832344279818,
 0.58984513209798939,
 1.0191773861094746,
 5.2424783744028893]